R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(data.table)
library(tidyr)
library(maps)
library(haven)
library(ggplot2)
library(dplyr)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

library(readxl)
hardship_complete <- read_excel("/Users/cristinacandido/Documents/Github/risk_wvs/code/Hardship_complete_2024.xlsx")
hardship_complete

NA
NA

hardship_complete$homiciderate=log(hardship_complete$homiciderate)
hardship_complete$gdp=log(hardship_complete$gdp)
hardship_complete$infantmortality=log(hardship_complete$infantmortality)
hardship_complete$lifeexpectancy=log(hardship_complete$lifeexpectancy)

hardship_complete

# Reverse Codierung
hardship_complete$homiciderate=scale(hardship_complete$homiciderate)
hardship_complete$gdp=scale(-hardship_complete$gdp)
hardship_complete$infantmortality=scale(hardship_complete$infantmortality)
hardship_complete$lifeexpectancy=scale(-hardship_complete$lifeexpectancy)
hardship_complete$gini=scale(hardship_complete$gini)
hardship_complete$femalemale_primedu=scale(-hardship_complete$femalemale_primedu)
hardship_complete

hardship_complete$hardship_index=(hardship_complete$homiciderate+hardship_complete$gdp+hardship_complete$gini+hardship_complete$lifeexpectancy+hardship_complete$infantmortality+hardship_complete$femalemale_primedu)/6

hardship_complete

# Laura added this code
library(dplyr)

# Impute missing values in the hardship indicators using the median
hardship_complete <- hardship_complete %>%
  mutate(across(c("homiciderate", "gdp", "infantmortality", "lifeexpectancy", "gini", "femalemale_primedu"),
                ~ifelse(is.na(.), median(., na.rm = TRUE), .)))

hardship_complete

# Use the mutate function to change the country name
hardship_complete <- hardship_complete %>%
  mutate(country = ifelse(label == "Serbia and Montenegro", "Serbia", label))
hardship_complete

# Data of Wave 5


WV5_data <- readRDS("/Users/cristinacandido/Documents/Github/risk_wvs/data/WVS/F00007944-WV5_Data_R_v20180912.rds")


# Convert WV5_data-object in data.frame 
WV5_data_df <- as.data.frame(WV5_data)

# show first five columns
WV5_data_df

#rename the variables
WV5_data <- WV5_data_df %>%
  rename(gender = V235, age = V237, country_code = V2, wave = V1, risktaking = V86, children = V56, married = V55, employed = V241, education = V238)
WV5_data

colnames(WV5_data)

  [1] "wave"          "V1A"           "V1B"           "country_code"  "V2A"           "V3"            "V4"            "V4_CO"         "V5"           
 [10] "V5_CO"         "V6"            "V6_CO"         "V7"            "V7_CO"         "V8"            "V8_CO"         "V9"            "V9_CO"        
 [19] "V10"           "V11"           "V12"           "V13"           "V14"           "V15"           "V16"           "V17"           "V18"          
 [28] "V19"           "V20"           "V21"           "V22"           "V23"           "V24"           "V25"           "V26"           "V27"          
 [37] "V28"           "V29"           "V30"           "V31"           "V32"           "V33"           "V34"           "V35"           "V36"          
 [46] "V37"           "V38"           "V39"           "V40"           "V41"           "V42"           "V43"           "V43_01"        "V43_02"       
 [55] "V43_03"        "V43_04"        "V43_05"        "V43_06"        "V43_07"        "V43_08"        "V43_09"        "V43_10"        "V43_11"       
 [64] "V43_12"        "V43_13"        "V43_14"        "V43_15"        "V43_16"        "V43_17"        "V43_18"        "V43_19"        "V43_20"       
 [73] "V43_21"        "V43_22"        "V43_23"        "V43_24"        "V43_25"        "V43_26"        "V43_27"        "V43_28"        "V43_29"       
 [82] "V43_30"        "V44"           "V45"           "V46"           "V47"           "V48"           "V49"           "V50"           "V51"          
 [91] "V52"           "V53"           "V54"           "married"       "children"      "V57"           "V58"           "V59"           "V60"          
[100] "V61"           "V62"           "V63"           "V64"           "V65"           "V66"           "V67"           "V68"           "V69"          
[109] "V69_HK"        "V70"           "V70_HK"        "V71"           "V72"           "V73"           "V73_HK"        "V74"           "V74_HK"       
[118] "V75"           "V76"           "V77"           "V78"           "V79"           "V80"           "V81"           "V82"           "V83"          
[127] "V84"           "V85"           "risktaking"    "V87"           "V88"           "V89"           "V90"           "V91"           "V92"          
[136] "V93"           "V94"           "V95"           "V96"           "V97"           "V98"           "V99"           "V100"          "V101"         
[145] "V102"          "V103"          "V104"          "V105"          "V106"          "V107"          "V108"          "V109"          "V110"         
[154] "V111"          "V112"          "V113"          "V114"          "V115"          "V116"          "V117"          "V118"          "V119"         
[163] "V120"          "V121"          "V122"          "V123"          "V124"          "V125"          "V126"          "V127"          "V128"         
[172] "V129"          "V130"          "V130_CA_1"     "V130_IQ_1"     "V130_IQ_2"     "V130_IQ_3"     "V130_IQ_4"     "V130_NZ_1"     "V130_NZ_2"    
[181] "V131"          "V132"          "V133"          "V134"          "V135"          "V136"          "V137"          "V138"          "V139"         
[190] "V140"          "V141"          "V142"          "V143"          "V144"          "V145"          "V146_00"       "V146_01"       "V146_02"      
[199] "V146_03"       "V146_04"       "V146_05"       "V146_06"       "V146_07"       "V146_08"       "V146_09"       "V146_10"       "V146_11"      
[208] "V146_12"       "V146_13"       "V146_14"       "V146_15"       "V146_16"       "V146_17"       "V146_18"       "V146_19"       "V146_20"      
[217] "V146_21"       "V146_22"       "V147"          "V148"          "V149"          "V150"          "V151"          "V151_IQ_A"     "V151_IQ_B"    
[226] "V152"          "V153"          "V154"          "V155"          "V156"          "V157"          "V158"          "V159"          "V160"         
[235] "V161"          "V162"          "V163"          "V164"          "V165"          "V166"          "V167"          "V168"          "V169"         
[244] "V170"          "V171"          "V172"          "V173"          "V174"          "V175"          "V176"          "V177"          "V178"         
[253] "V179"          "V180"          "V181"          "V182"          "V183"          "V184"          "V185"          "V186"          "V187"         
[262] "V188"          "V189"          "V190"          "V191"          "V192"          "V193"          "V194"          "V195"          "V196"         
[271] "V197"          "V198"          "V199"          "V200"          "V201"          "V202"          "V203"          "V204"          "V205"         
[280] "V206"          "V207"          "V208"          "V209"          "V210"          "V211"          "V212"          "V213A"         "V213B"        
[289] "V213C"         "V213D"         "V213E"         "V213F"         "V213G"         "V213H"         "V213K"         "V213L"         "V213M"        
[298] "V213N"         "V214"          "V215"          "V216"          "V217"          "V218"          "V219"          "V220"          "V221"         
[307] "V222"          "V223"          "V224"          "V225"          "V226"          "V227"          "V228"          "V229"          "V230"         
[316] "V231"          "V232"          "V233"          "V233A"         "V234"          "gender"        "V236"          "age"           "education"    
[325] "V238CS"        "V239"          "V240"          "employed"      "V242"          "V242A_CO"      "V243"          "V244"          "V245"         
[334] "V246"          "V247"          "V248"          "V249"          "V250"          "V251"          "V252"          "V252B"         "V253"         
[343] "V253CS"        "V254"          "V255"          "V255CS"        "V256"          "V257"          "V257B"         "V257C"         "V258"         
[352] "V259"          "V259A"         "V260"          "V261"          "V262"          "V263"          "V264"          "V265"          "S024"         
[361] "S025"          "Y001"          "Y002"          "Y003"          "SACSECVAL"     "SECVALWGT"     "RESEMAVAL"     "WEIGHTB"       "I_AUTHORITY"  
[370] "I_NATIONALISM" "I_DEVOUT"      "DEFIANCE"      "WEIGHT1A"      "I_RELIGIMP"    "I_RELIGBEL"    "I_RELIGPRAC"   "DISBELIEF"     "WEIGHT2A"     
[379] "I_NORM1"       "I_NORM2"       "I_NORM3"       "RELATIVISM"    "WEIGHT3A"      "I_TRUSTARMY"   "I_TRUSTPOLICE" "I_TRUSTCOURTS" "SCEPTICISM"   
[388] "WEIGHT4A"      "I_INDEP"       "I_IMAGIN"      "I_NONOBED"     "AUTONOMY"      "WEIGHT1B"      "I_WOMJOB"      "I_WOMPOL"      "I_WOMEDU"     
[397] "EQUALITY"      "WEIGHT2B"      "I_HOMOLIB"     "I_ABORTLIB"    "I_DIVORLIB"    "CHOICE"        "WEIGHT3B"      "I_VOICE1"      "I_VOICE2"     
[406] "I_VOI2_00"     "VOICE"         "WEIGHT4B"      "S001"          "S007"          "S018"          "S019"          "S021"          "COW"          

#select only the variables of interest
WV5_data <- WV5_data %>%
  dplyr::select(gender, age, country_code, wave, risktaking, children, married, employed, education)
WV5_data

#Read Dataset (Wave 6)

WV6_data <- load("/Users/cristinacandido/Documents/Github/risk_wvs/data/WVS/WV6_Data_R_v20201117.rdata") 
WV6_data <- WV6_Data_R_v20201117 
print(WV6_data)

#Read Dataset (Wave 6)

WV6_data <- load("/Users/cristinacandido/Documents/Github/risk_wvs/data/WVS/WV6_Data_R_v20201117.rdata") 
WV6_data <- WV6_Data_R_v20201117 
print(WV6_data)

WV6_data <- WV6_data %>%
  rename(wave = V1, gender = V240, age = V242,country_code = V2, risktaking = V76, children = V58, married = V57, employed = V229, education = V248)


#select only the variables of interest

WV6_data <- WV6_data %>%
  dplyr::select(wave, gender, age, country_code,risktaking, children, married, employed, education)
WV6_data

countrynames = read.csv("/Users/cristinacandido/Documents/Github/risk_wvs/data/WVS/countrynames.txt", header=FALSE,as.is=TRUE)
colnames(countrynames) = c("code", "name")
WV6_data$country = countrynames$name [match(WV6_data$country_code, countrynames$code)]
table(WV6_data$country)
WV6_data

WVS_data = rbind(WV5_data, WV6_data)
WVS_data

country_counts <- WVS_data %>%
  count(country)

# Print the result
print(country_counts)

NA
NA
NA
NA
NA
NA
NA
NA
NA
NA

WVS_data = subset(WVS_data, risktaking > 0 & gender > 0 & age >0 & education > 0 & employed > 0 & married > 0 & children >= 0)
data_Wave5 = subset(WV5_data, risktaking > 0 & gender > 0 & age >0 & education > 0 & employed > 0 & married > 0 & children >= 0)
data_Wave6 = subset(WV6_data, risktaking > 0 & gender > 0 & age >0 & education > 0 & employed > 0 & married > 0 & children >= 0)
WVS_data <- na.omit(WVS_data)
data_Wave5 <- na.omit(data_Wave5)
data_Wave6 <- na.omit(data_Wave6)









# Use the mutate function to change the country name
WVS_data <- WVS_data %>%
  mutate(country = ifelse(country == "Great Britain", "United Kingdom", country))

unique(WVS_data$country)

 [1] "Andorra"             "Argentina"           "Australia"           "Brazil"              "Bulgaria"           
 [6] "Burkina Faso"        "Canada"              "Cyprus (G)"          "Chile"               "China"              
[11] "Egypt"               "Ethiopia"            "Finland"             "France"              "Georgia"            
[16] "Germany"             "Ghana"               "Hungary"             "India"               "Indonesia"          
[21] "Iran"                "Japan"               "Malaysia"            "Mali"                "Mexico"             
[26] "Moldova"             "Morocco"             "Netherlands"         "Norway"              "Peru"               
[31] "Poland"              "Romania"             "Russia"              "Rwanda"              "Slovenia"           
[36] "South Africa"        "South Korea"         "Spain"               "Sweden"              "Switzerland"        
[41] "Taiwan"              "Thailand"            "Trinidad and Tobago" "Turkey"              "Ukraine"            
[46] "United Kingdom"      "Uruguay"             "Viet Nam"            "Zambia"              "Algeria"            
[51] "Armenia"             "Azerbaijan"          "Belarus"             "Colombia"            "Ecuador"            
[56] "Estonia"             "Haiti"               "Iraq"                "Jordan"              "Kazakhstan"         
[61] "Kuwait"              "Kyrgyzstan"          "Lebanon"             "Libya"               "New Zealand"        
[66] "Nigeria"             "Pakistan"            "Palestine"           "Philippines"         "Qatar"              
[71] "Singapore"           "Tunisia"             "United States"       "Uzbekistan"          "Yemen"              
[76] "Zimbabwe"           

# age, sex 
model1 = lmer(T_score_risktaking ~ 1 +scale(age)+factor(gender) + (1+scale(age)+factor(gender)|country),data = WVS_mixed_model, control = lmerControl(optimizer = "bobyqa"))
summary_model1=summary(model1)
summary_model1

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: T_score_risktaking ~ 1 + scale(age) + factor(gender) + (1 + scale(age) +      factor(gender) | country)
   Data: WVS_mixed_model
Control: lmerControl(optimizer = "bobyqa")

REML criterion at convergence: 1081364

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4737 -0.7821 -0.0770  0.7536  3.2222 

Random effects:
 Groups   Name            Variance Std.Dev. Corr     
 country  (Intercept)      6.6199  2.5729            
          scale(age)       0.7621  0.8730   0.30     
          factor(gender)1  0.8747  0.9353   0.13 0.27
 Residual                 84.6046  9.1981            
Number of obs: 148527, groups:  country, 76

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      51.4031     0.2977 75.2144  172.66   <2e-16 ***
scale(age)       -2.0093     0.1042 73.1211  -19.28   <2e-16 ***
factor(gender)1  -2.3174     0.1194 72.3693  -19.41   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(g)
scale(age)  0.291        
fctr(gndr)1 0.072  0.236 

#model 2
library(lme4)

library(lme4)
library(lmerTest)

# Define the lmer model and assign it to 'model_2'
model_2 <- lmer(T_score_risktaking ~ 1 + scale(age) + factor(gender) + factor(children) + 
                 factor(married) + factor(education) + factor(employed) +
                 (1 + scale(age) + factor(gender) + factor(children) + factor(employed) +
                  factor(married) | country),
                data = WVS_mixed_model, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))


# Display the summary of the model
summary_model_2 = summary(model_2)
summary_model_2

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: T_score_risktaking ~ 1 + scale(age) + factor(gender) + factor(children) +  
    factor(married) + factor(education) + factor(employed) +  
    (1 + scale(age) + factor(gender) + factor(children) + factor(employed) +          factor(married) | country)
   Data: WVS_mixed_model
Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

REML criterion at convergence: 1079996

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5457 -0.7815 -0.0813  0.7425  3.1678 

Random effects:
 Groups   Name              Variance Std.Dev. Corr                         
 country  (Intercept)        5.4671  2.3382                                
          scale(age)         0.5367  0.7326    0.24                        
          factor(gender)1    0.9298  0.9642    0.03  0.24                  
          factor(children)1  0.7895  0.8885    0.07  0.19  0.11            
          factor(employed)1  0.3133  0.5597    0.04  0.03  0.02 -0.27      
          factor(married)1   0.4120  0.6419    0.22  0.43  0.57  0.20 -0.13
 Residual                   83.7217  9.1500                                
Number of obs: 148527, groups:  country, 76

Fixed effects:
                     Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)         5.204e+01  2.809e-01  8.000e+01 185.238  < 2e-16 ***
scale(age)         -1.478e+00  9.081e-02  7.228e+01 -16.275  < 2e-16 ***
factor(gender)1    -2.171e+00  1.234e-01  7.259e+01 -17.597  < 2e-16 ***
factor(children)1  -1.274e+00  1.284e-01  7.199e+01  -9.920 4.14e-15 ***
factor(married)1   -8.367e-01  9.872e-02  6.267e+01  -8.475 5.49e-12 ***
factor(education)1  7.809e-01  6.489e-02  1.407e+05  12.033  < 2e-16 ***
factor(employed)1   1.056e-01  8.590e-02  7.071e+01   1.229    0.223    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(g) fctr(g)1 fctr(c)1 fctr(mr)1 fctr(d)1
scale(age)   0.224                                            
fctr(gndr)1 -0.022  0.212                                     
fctr(chld)1 -0.016  0.047  0.032                              
fctr(mrrd)1  0.134  0.297  0.409   -0.083                     
fctr(dctn)1 -0.172  0.081  0.013    0.026   -0.013            
fctr(mply)1 -0.029  0.061  0.083   -0.203   -0.099    -0.063  

WVS_mixed_model

NA
NA

# age, sex 
model1 = lmer(T_score_risktaking ~ 1 +scale(age)+factor(gender) + (1+scale(age)+factor(gender)|country),data = WVS_mixed_model, control = lmerControl(optimizer = "bobyqa"))
summary_model1=summary(model1)
summary_model1

model_3 <- lmer(T_score_risktaking ~ 1 + scale(age) * hardship_index + 
                    factor(gender) * hardship_index + factor(married) + factor(children) + 
                    factor(education) + factor(employed) + 
                    (1 + scale(age) + factor(married) + factor(children) + 
                     factor(education) + factor(employed) | country),
                data = WVS_mixed_model,control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

summary_model_3 = summary(model_3) 

summary_model_3

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: T_score_risktaking ~ 1 + scale(age) * hardship_index + factor(gender) *  
    hardship_index + factor(married) + factor(children) + factor(education) +  
    factor(employed) + (1 + scale(age) + factor(married) + factor(children) +      factor(education) + factor(employed) | country)
   Data: WVS_mixed_model
Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

REML criterion at convergence: 1080093

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.57571 -0.78110 -0.08036  0.74750  3.13409 

Random effects:
 Groups   Name               Variance Std.Dev. Corr                         
 country  (Intercept)         5.5510  2.3561                                
          scale(age)          0.3952  0.6286    0.17                        
          factor(married)1    0.4028  0.6347    0.13  0.13                  
          factor(children)1   0.8063  0.8980    0.09  0.06  0.25            
          factor(education)1  0.6332  0.7957   -0.20  0.03  0.07 -0.13      
          factor(employed)1   0.2881  0.5368    0.00  0.08 -0.29 -0.25  0.16
 Residual                    83.8078  9.1547                                
Number of obs: 148527, groups:  country, 76

Fixed effects:
                                 Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                     5.204e+01  2.862e-01  7.433e+01 181.838  < 2e-16 ***
scale(age)                     -1.442e+00  8.035e-02  7.109e+01 -17.949  < 2e-16 ***
hardship_index                  4.588e-01  3.840e-01  7.514e+01   1.195    0.236    
factor(gender)1                -2.081e+00  5.032e-02  1.285e+05 -41.359  < 2e-16 ***
factor(married)1               -8.145e-01  9.857e-02  6.114e+01  -8.263 1.53e-11 ***
factor(children)1              -1.302e+00  1.292e-01  7.512e+01 -10.070 1.39e-15 ***
factor(education)1              7.768e-01  1.205e-01  5.044e+01   6.447 4.24e-08 ***
factor(employed)1               1.031e-01  8.357e-02  6.975e+01   1.233    0.222    
scale(age):hardship_index       4.751e-01  1.122e-01  7.216e+01   4.233 6.68e-05 ***
hardship_index:factor(gender)1  5.799e-01  6.902e-02  6.340e+04   8.403  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(g) hrdsh_ fctr(g)1 fctr(mr)1 fctr(c)1 fctr(d)1 fctr(mp)1 sc():_
scale(age)   0.158                                                                    
hardshp_ndx  0.067  0.034                                                             
fctr(gndr)1 -0.098  0.042  0.002                                                      
fctr(mrrd)1  0.066  0.078 -0.005  0.044                                               
fctr(chld)1  0.005 -0.062 -0.003 -0.087   -0.054                                      
fctr(dctn)1 -0.285  0.080  0.014  0.016    0.030    -0.069                            
fctr(mply)1 -0.056  0.095  0.006  0.145   -0.190    -0.186    0.046                   
scl(g):hrd_  0.032  0.091  0.191  0.002   -0.005     0.005   -0.002   -0.035          
hrdshp_:()1 -0.002  0.011 -0.092  0.057   -0.012    -0.006    0.006    0.030     0.023

random_effects <- ranef(model_2)
random_effects

$country
                    (Intercept)   scale(age) factor(gender)1 factor(children)1 factor(employed)1 factor(married)1
Algeria              1.93351855 -0.854619155      0.87619718      -2.437881334        0.17698875     -0.377226443
Andorra              1.91309921  0.100909057     -0.66981587      -0.632350380        0.20524734     -0.686912835
Argentina           -1.57325789  0.057299092     -0.25780065      -1.067941885       -0.04621841     -0.553872390
Armenia              1.08284389 -0.454274344     -0.89771357      -0.834967062        0.17426467     -0.250770715
Australia           -0.14386428 -0.471375287     -1.06538548      -0.103260090        0.34608388     -0.365152174
Azerbaijan           0.45573171  0.004004243     -2.44218293       0.082693851       -0.21915868     -0.476480146
Belarus             -0.11536490 -1.000016058     -1.08120720      -0.467055853       -0.59713729     -0.046875401
Brazil              -2.84287325  0.262612027     -0.36729059      -0.148830175       -0.34944728     -0.655399532
Bulgaria            -0.77688584 -0.074197571      0.35410308       0.108171529        0.66545761     -0.168532781
Burkina Faso         1.61422833  1.156372503     -0.15049497       0.527033659       -0.25535458      0.423175478
Canada               0.39907784 -0.397232782     -1.02665629      -0.204786757        0.83298979     -0.733722779
Chile                1.03350248 -0.644167015      0.06242265      -1.191152814        0.30273676      0.182724595
China               -2.09375223 -0.445731013      0.50663419      -0.603547862       -0.41625835     -0.223192477
Colombia             0.22689773  0.422663335     -0.91548157      -0.850065668        0.05700229     -0.046926087
Cyprus (G)           2.25496876 -0.776460854     -0.79745197       0.050374871        0.20526698     -0.994122261
Ecuador              2.25294069 -0.040563584      0.75768462       1.082273600       -0.18216772      0.606966414
Egypt               -3.93763995  0.877857068      1.47333952      -0.977664255        0.69617277      0.506927589
Estonia             -1.06399862 -0.805260418     -0.20652339      -0.596528530        0.04115956     -0.263441563
Ethiopia             0.32293567  0.946455342      1.98005066      -0.482906872        1.19465275      0.520837889
Finland             -2.30711075 -0.193362752      0.06835226       0.109995921       -0.09624712     -0.446129133
France              -0.03262671 -0.315869185     -0.97254579      -0.756893462       -0.18841330     -0.721010031
Georgia             -0.34804857  0.313082762     -1.59932353       0.785006390        0.24885966     -0.322662614
Germany             -2.81688356 -0.866667483     -0.56600732      -0.067064088        0.10227997     -0.540483663
Ghana                4.20166977  0.817382387      1.03493510       1.320604468        0.06952553      0.935149120
Haiti               -6.53977492 -1.391127796      1.91289649       1.468940014       -0.89280260      0.533482158
Hungary             -1.99653091 -0.996752173     -0.17049773      -0.160714296        0.07124129     -0.186539589
India                1.43482568  0.617405603      0.45019893       1.836342889       -1.09159157      0.408036403
Indonesia            2.91129285  0.717729512     -0.37251725       0.036310599        0.61487372      0.351255947
Iran                 0.25704414 -0.273465758      0.84529203      -0.730085449        0.58356967     -0.185573571
Iraq                -1.06492016  0.012316662      0.43906584       0.138047194       -0.32549734      0.649008088
Japan               -6.33127212  0.980885912     -0.11206958       0.564989629       -0.44804909     -0.169376931
Jordan               1.80284181 -0.238660302     -0.04733756      -0.275114119        0.13106070      0.079727138
Kazakhstan          -2.61497176  0.099878889      0.38099799      -0.441649153        0.38909549      0.080016173
Kuwait               2.09145419 -1.486900320      0.10785659       0.216244412       -0.71982856     -0.372559647
Kyrgyzstan           0.61067751  1.072544199      1.07198178       0.701056752       -0.33081703      0.900953951
Lebanon              2.84470796 -0.023345343      0.08714121       0.523974067       -0.06162559      0.535210026
Libya                0.80327073 -0.475177836     -0.57060729       0.015373078        0.24676530     -0.259391475
Malaysia            -2.02257329  0.667498135      0.69083376       0.230243197       -0.12134511      0.648213460
Mali                 2.81563100  1.078784421      1.36765639       0.245910074       -0.30999444      1.287801120
Mexico              -0.40013922  0.561817933     -0.55226850      -1.131405064        0.01974498     -0.487453888
Moldova             -1.96423291 -0.329843137      1.07398956      -0.873396236        0.07687739      0.075323668
Morocco              0.68539027  0.249483178     -0.44534432      -0.309370915        0.01082943     -0.178509821
Netherlands         -2.06588418 -0.228085056     -0.89179246       0.412747451       -0.17127304     -0.633635212
New Zealand          0.60732681 -0.774713993      0.04302897      -0.175577939       -0.11512608     -0.279411963
Nigeria              5.59377567  1.275958637      1.62013094       0.389674844        0.37813269      1.361593143
Norway               0.99797819 -0.771979285      0.88446601      -0.045167601       -0.33921268     -0.013483356
Pakistan             4.22659980  1.325790041      0.51442825       0.928501617       -0.97792666      0.666654501
Palestine            0.03194219 -0.494835906     -0.14515969       0.292953341        0.50604479     -0.015691232
Peru                -1.61516779  0.633332407     -0.68009631      -0.717362556       -0.39962778     -0.181622730
Philippines          3.70377389  0.943369150      1.13088110       0.747509770       -0.08964864      0.356984425
Poland               2.19905667  0.037084593     -0.45404816      -0.900087529        0.25558657      0.232224025
Qatar                0.86028228  0.502376889     -0.94153784       0.200394431        0.45200189     -0.086160547
Romania             -2.52089552 -0.155453441     -0.24318156      -0.213424238        0.59038608     -0.046222935
Russia              -0.03564662  0.201433455     -0.04277014       0.338525048        0.55241645     -0.367819701
Rwanda               0.83617515  1.171444046      1.19476510       1.697279263       -0.29018529      0.335814980
Singapore            2.17597132  0.591444549      0.38684581      -0.017294940       -0.08165496      0.363947899
Slovenia             0.21102243 -1.245719908     -0.62662173      -0.617380784        0.28218849     -1.002740221
South Africa         2.37631040  0.311956458      1.26337006       1.080146122        0.68094981      0.751652547
South Korea          1.36893323 -0.405664050      0.80080531       0.341944092       -0.04931702      0.322300385
Spain                0.82809176 -0.466229484      0.62854866      -0.304041861        0.09422269     -0.358645793
Sweden               0.09920433 -0.328664254     -0.47169418      -0.328195601       -0.42679317     -0.438405460
Switzerland         -1.56290696  0.084346338     -1.48036381       0.605515233        0.14178415     -0.884764664
Taiwan              -4.35319381 -0.569177868     -0.09471962      -0.003120567       -0.38488863     -0.220827542
Thailand            -0.09358157  0.844688730      1.30128241       0.253396867       -0.53318528      1.057267847
Trinidad and Tobago  0.24160167 -0.343056024     -0.36770149       1.121472297        0.13881740      0.059010394
Tunisia             -0.42814434 -1.362242692     -0.53278180      -0.600878765       -0.28277064     -0.244015679
Turkey               0.94251912  0.120341608      0.37598170      -0.387346063        0.29373891      0.028586139
Ukraine             -1.58175919 -0.668879956     -0.36649799      -0.519892277        0.26972976     -0.001257245
United Kingdom       1.15668394 -0.030727372     -1.56593196      -0.489635482        0.10560093     -0.751106765
United States       -0.70200525 -0.752298720     -0.40433841       0.416924321        0.25359660     -0.253434481
Uruguay             -3.20629770 -0.108243534      0.34971574      -0.439094936       -0.31919572     -0.167863204
Uzbekistan           1.00109484  0.687012035     -0.28533330       0.716274322        0.20755393      0.141926144
Viet Nam            -3.44301278  0.598580704     -0.35309866       0.622440999       -0.56204886      0.334138408
Yemen               -3.91719513  0.062896576     -0.10339324      -0.646218561       -0.08885691     -0.219935374
Zambia               1.41695079  0.455513249      0.41766478       1.043367849       -0.57988565      1.261564524
Zimbabwe             1.68856739  0.396489975     -1.11595904       0.496697944       -0.32194635     -0.119112532

with conditional variances for “country” 

# Extract random effects for 'country'
# Assuming 'model' is your fitted lmer model
random_effects <- ranef(model_2)
random_effects

colnames(random_effects)

NULL

coefsallmodels=rbind(summary_model1$coefficients,
summary_model_2$coefficients,
summary_model_3$coefficients[c(1:2,4:8,3,9:10),])

write.csv(coefsallmodels,"coefsallmodels.csv")

# Read the CSV file into a data frame
gbd_mentalhealth <- read_excel("/Users/cristinacandido/Documents/Github/risk_wvs/GBD_mentalhealth.xlsx")

gbd_mentalhealth

#select only the variables of interest
gbd_mentalhealth <- gbd_mentalhealth %>%
  dplyr::select(country, gender, age, cause, val, Measure)
gbd_mentalhealth


library(dplyr)

# Group data by country and age group, and calculate summary statistics
summary_by_country_age <- gbd_mentalhealth %>%
  group_by(country, age) %>%
  summarise(
    mean_DALYs = mean(val))  # Calculate mean of DALYs

library(dplyr)

# Assuming 'summary_by_country_age' contains your summarized dataset
mean_by_country <- summary_by_country_age %>%
  group_by(country) %>%
  summarise(mean_DALYs = mean(mean_DALYs))

# View the resulting mean by country
print(mean_by_country)

#log transform
mean_by_country$mean_DALYs=log(mean_by_country$mean_DALYs)
 
mean_by_country 

#Reverse codierung 
mean_by_country$mean_DALYs=scale(mean_by_country$mean_DALYs)

mean_by_country

#rename mean_DALYS
mental_health_index <- mean_by_country %>%
  rename('mental_health' = mean_DALYs)
mental_health_index

library(dplyr)

#######Anxiety disorders#########
#Filter data for a specific mental disorder (e.g., Anxiety Disorders)
anxiety_disorders <- gbd_mentalhealth %>%
  filter(cause == "Anxiety disorders")  # Change "Anxiety Disorders" to the desired disorder

# Calculate the mean of 'val' (Disability-Adjusted Life Years) across locations and ages
mean_DALYs <- mean(anxiety_disorders$val)

# Optionally, if you want to calculate mean by country:
anxiety_disorders <- anxiety_disorders %>%
  group_by(country) %>%
  summarise(mean_Anxiety_disorders = mean(val))

anxiety_disorders

#log transform
anxiety_disorders$mean_Anxiety_disorders=log(anxiety_disorders$mean_Anxiety_disorders)
 
anxiety_disorders

#Reverse codierung 
anxiety_disorders$mean_Anxiety_disorders=scale(anxiety_disorders$mean_Anxiety_disorders)
anxiety_disorders
mental_health_index

gbd_mentalhealth

# Assuming gbd_mentalhealth is your data frame containing the 'cause' column

# Get unique values of 'cause' column
unique_causes <- unique(gbd_mentalhealth$cause)

# Create a data frame with unique causes
cause_table <- data.frame(Cause = unique_causes)

# Print the cause table
print(cause_table)

#######Bulimia nervosa########
library(dplyr)

# Filter the data for 'Bulimia Nervosa'
bulimia_nervosa <- gbd_mentalhealth %>%
  filter(cause == "Bulimia nervosa")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(bulimia_nervosa$val)

# Optionally, calculate mean by country
bulimia_nervosa <- bulimia_nervosa %>%
  group_by(country) %>%
  summarise(mean_bulimia_nervosa = mean(val))

# Log transform the mean_bulimia_nervosa column
bulimia_nervosa$mean_bulimia_nervosa <- log(bulimia_nervosa$mean_bulimia_nervosa)

# Scale the mean_bulimia_nervosa column (if needed)
bulimia_nervosa$mean_bulimia_nervosa <- scale(bulimia_nervosa$mean_bulimia_nervosa)

# Print the resulting data frame
print(bulimia_nervosa)
print(anxiety_disorders)

gbd_mentalhealth

#####Attention-deficit/hyperactivity disorder
library(dplyr)


ADHD <- gbd_mentalhealth %>%
  filter(cause == "Attention-deficit/hyperactivity disorder")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(ADHD$val)

# Optionally, calculate mean by country
ADHD <- ADHD %>%
  group_by(country) %>%
  summarise(mean_ADHD = mean(val))

# Log transform the mean_bulimia_nervosa column
ADHD$mean_ADHD <- log(ADHD$mean_ADHD)

# Scale the mean_bulimia_nervosa column (if needed)
ADHD$mean_ADHD <- scale(ADHD$mean_ADHD)

# Print the resulting data frame
print(ADHD)

#########Idiopathic development intellectual ability 
library(dplyr)


Idiopathic_developmental_intellectual_disability <- gbd_mentalhealth %>%
  filter(cause == "Idiopathic developmental intellectual disability")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(Idiopathic_developmental_intellectual_disability$val)

# Optionally, calculate mean by country
Idiopathic_developmental_intellectual_disability <- Idiopathic_developmental_intellectual_disability %>%
  group_by(country) %>%
  summarise(mean_Idiopathic_developmental_intellectual_disability = mean(val))

# Log transform the mean_bulimia_nervosa column
Idiopathic_developmental_intellectual_disability$mean_Idiopathic_developmental_intellectual_disability <- log(Idiopathic_developmental_intellectual_disability$mean_Idiopathic_developmental_intellectual_disability)

# Scale the mean_bulimia_nervosa column (if needed)
Idiopathic_developmental_intellectual_disability$mean_Idiopathic_developmental_intellectual_disability <- scale(Idiopathic_developmental_intellectual_disability$mean_Idiopathic_developmental_intellectual_disability)

# Print the resulting data frame
print(Idiopathic_developmental_intellectual_disability)

######Anorexia Nervosa

anorexia_nervosa <- gbd_mentalhealth %>%
  filter(cause == "Anorexia nervosa")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(anorexia_nervosa$val)

# Optionally, calculate mean by country
anorexia_nervosa <- anorexia_nervosa %>%
  group_by(country) %>%
  summarise(mean_anorexia_nervosa = mean(val))

# Log transform the mean_bulimia_nervosa column
anorexia_nervosa$mean_anorexia_nervosa <- log(anorexia_nervosa$mean_anorexia_nervosa)

# Scale the mean_bulimia_nervosa column (if needed)
anorexia_nervosa$mean_anorexia_nervosa <- scale(anorexia_nervosa$mean_anorexia_nervosa)

# Print the resulting data frame
print(anorexia_nervosa)

#####Depressive disorders#######

depressive_disorders <- gbd_mentalhealth %>%
  filter(cause == "Depressive disorders")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(depressive_disorders$val)

# Optionally, calculate mean by country
depressive_disorders <- depressive_disorders %>%
  group_by(country) %>%
  summarise(mean_depressive_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
depressive_disorders$mean_depressive_disorders <- log(depressive_disorders$mean_depressive_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
depressive_disorders$mean_depressive_disorders <- scale(depressive_disorders$mean_depressive_disorders)

# Print the resulting data frame
print(depressive_disorders)

#######Autismus spectrum disorders######

autismus_spectrum_disorders <- gbd_mentalhealth %>%
  filter(cause == "Autism spectrum disorders")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(autismus_spectrum_disorders$val)

# Optionally, calculate mean by country
autismus_spectrum_disorders <- autismus_spectrum_disorders %>%
  group_by(country) %>%
  summarise(mean_autismus_spectrum_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
autismus_spectrum_disorders$mean_autismus_spectrum_disorders <- log(autismus_spectrum_disorders$mean_autismus_spectrum_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
autismus_spectrum_disorders$mean_autismus_spectrum_disorders <- scale(autismus_spectrum_disorders$mean_autismus_spectrum_disorders)

# Print the resulting data frame
print(autismus_spectrum_disorders)

######Schizophrenia#######
schizophrenia <- gbd_mentalhealth %>%
  filter(cause == "Schizophrenia")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(schizophrenia$val)

# Optionally, calculate mean by country
schizophrenia <- schizophrenia %>%
  group_by(country) %>%
  summarise(mean_schizophrenia = mean(val))

# Log transform the mean_bulimia_nervosa column
schizophrenia$mean_schizophrenia <- log(schizophrenia$mean_schizophrenia)

# Scale the mean_bulimia_nervosa column (if needed)
schizophrenia$mean_schizophrenia <- scale(schizophrenia$mean_schizophrenia)

# Print the resulting data frame
print(schizophrenia)

#######Conduct disorders#########
conduct_disorders <- gbd_mentalhealth %>%
  filter(cause == "Conduct disorder")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(conduct_disorders$val)

# Optionally, calculate mean by country
conduct_disorders <- conduct_disorders %>%
  group_by(country) %>%
  summarise(mean_conduct_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
conduct_disorders$mean_conduct_disorders <- log(conduct_disorders$mean_conduct_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
conduct_disorders$mean_conduct_disorders <- scale(conduct_disorders$mean_conduct_disorders)

# Print the resulting data frame
print(conduct_disorders)

########Eating disorders#########
eating_disorders <- gbd_mentalhealth %>%
  filter(cause == "Eating disorders")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(eating_disorders$val)

# Optionally, calculate mean by country
eating_disorders <- eating_disorders %>%
  group_by(country) %>%
  summarise(mean_eating_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
eating_disorders$mean_eating_disorders <- log(eating_disorders$mean_eating_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
eating_disorders$mean_eating_disorders <- scale(eating_disorders$mean_eating_disorders)

# Print the resulting data frame
print(eating_disorders)

########Bipolar disorder#########
bipolar_disorder <- gbd_mentalhealth %>%
  filter(cause == "Bipolar disorder")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(bipolar_disorder$val)

# Optionally, calculate mean by country
bipolar_disorder <- bipolar_disorder %>%
  group_by(country) %>%
  summarise(mean_bipolar_disorder = mean(val))

# Log transform the mean_bulimia_nervosa column
bipolar_disorder$mean_bipolar_disorder <- log(bipolar_disorder$mean_bipolar_disorder)

# Scale the mean_bulimia_nervosa column (if needed)
bipolar_disorder$mean_bipolar_disorder <- scale(bipolar_disorder$mean_bipolar_disorder)

# Print the resulting data frame
print(bipolar_disorder)

print(cause_table)

########Substance use disorders########

substance_use_disorders <- gbd_mentalhealth %>%
  filter(cause == "Substance use disorders")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(substance_use_disorders$val)

# Optionally, calculate mean by country
substance_use_disorders <- substance_use_disorders %>%
  group_by(country) %>%
  summarise(mean_substance_use_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
substance_use_disorders$mean_substance_use_disorders <- log(substance_use_disorders$mean_substance_use_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
substance_use_disorders$mean_substance_use_disorders <- scale(substance_use_disorders$mean_substance_use_disorders)

# Print the resulting data frame
print(substance_use_disorders)

####Drug use disorders

drug_use_disorders <- gbd_mentalhealth %>%
  filter(cause == "Drug use disorders")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(drug_use_disorders$val)

# Optionally, calculate mean by country
drug_use_disorders <- drug_use_disorders %>%
  group_by(country) %>%
  summarise(mean_drug_use_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
drug_use_disorders$mean_drug_use_disorders <- log(drug_use_disorders$mean_drug_use_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
drug_use_disorders$mean_drug_use_disorders <- scale(drug_use_disorders$mean_drug_use_disorders)

# Print the resulting data frame
print(drug_use_disorders)

##########Alcohol use disorders
alcohol_use_disorders <- gbd_mentalhealth %>%
  filter(cause == "Alcohol use disorders")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(alcohol_use_disorders$val)

# Optionally, calculate mean by country
alcohol_use_disorders <- alcohol_use_disorders %>%
  group_by(country) %>%
  summarise(mean_alcohol_use_disorders = mean(val))

# Log transform the mean_bulimia_nervosa column
alcohol_use_disorders$mean_alcohol_use_disorders <- log(alcohol_use_disorders$mean_alcohol_use_disorders)

# Scale the mean_bulimia_nervosa column (if needed)
alcohol_use_disorders$mean_alcohol_use_disorders <- scale(alcohol_use_disorders$mean_alcohol_use_disorders)

# Print the resulting data frame
print(alcohol_use_disorders)

#######Major depressive disorders
major_depressive_disorder <- gbd_mentalhealth %>%
  filter(cause == "Major depressive disorder")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(major_depressive_disorder$val)

# Optionally, calculate mean by country
major_depressive_disorder <- major_depressive_disorder %>%
  group_by(country) %>%
  summarise(mean_major_depressive_disorder = mean(val))

# Log transform the mean_bulimia_nervosa column
major_depressive_disorder$mean_major_depressive_disorder <- log(major_depressive_disorder$mean_major_depressive_disorder)

# Scale the mean_bulimia_nervosa column (if needed)
major_depressive_disorder$mean_major_depressive_disorder <- scale(major_depressive_disorder$mean_major_depressive_disorder)

# Print the resulting data frame
print(major_depressive_disorder)

#######Dysthymia########

dysthymia <- gbd_mentalhealth %>%
  filter(cause == "Dysthymia")

# Calculate the mean of 'val' (DALYs) across all locations
mean_DALYs <- mean(dysthymia$val)

# Optionally, calculate mean by country
dysthymia <- dysthymia %>%
  group_by(country) %>%
  summarise(mean_dysthymia = mean(val))

# Log transform the mean_bulimia_nervosa column
dysthymia$mean_dysthymia <- log(dysthymia$mean_dysthymia)

# Scale the mean_bulimia_nervosa column (if needed)
dysthymia$mean_dysthymia <- scale(dysthymia$mean_dysthymia)

# Print the resulting data frame
print(dysthymia)

library(dplyr)

# Perform left joins for each data frame on 'country'
final_mental_health_index <- mental_health_index %>%
  left_join(bulimia_nervosa, by = "country") %>%
  left_join(ADHD, by = "country") %>%
  left_join(Idiopathic_developmental_intellectual_disability, by = "country") %>%
  left_join(anorexia_nervosa, by = "country") %>%
  left_join(depressive_disorders, by = "country") %>%
  left_join(autismus_spectrum_disorders, by = "country") %>%
  left_join(conduct_disorders, by = "country") %>%
  left_join(schizophrenia, by = "country") %>%
  left_join(eating_disorders, by = "country") %>%
  left_join(bipolar_disorder, by = "country") %>%
  left_join(drug_use_disorders, by = "country") %>%
  left_join(alcohol_use_disorders, by = "country") %>%
  left_join(substance_use_disorders, by = "country") %>%
  left_join(major_depressive_disorder, by = "country") %>%
  left_join(dysthymia, by = "country")

# Display the resulting data frame
print(final_mental_health_index)

# Show the first few rows of the resulting data frame
head(final_mental_health_index)

indicators <- left_join(final_mental_health_index, hardship, by = "country")
indicators
head(indicators)

new_data <- left_join (WVS_data, indicators, by = "country")
new_data


# Transfrom risk item such that high values represent more risk taking
new_data$risktaking = 6 - new_data$risktaking + 1

  
# Transform risk variable into T-score (mean = 50, sd = 10)
new_data$T_score_risktaking = 10*scale(new_data$risktaking, center=TRUE,scale=TRUE)+50

new_data

#Transform risk variable into Z score 

# Assuming T-scores have a mean of 50 and a standard deviation of 10
new_data$Z_score_risktaking = (new_data$T_score_risktaking - 50) / 10

# Print the resulting data frame
print(new_data)

new_data <- new_data %>%
  group_by(country) %>%
  mutate(z_score_age = scale(age))
new_data

library(lme4)

library(lme4)

model <- lmer(T_score_risktaking ~ scale(z_score_age) * mental_health +
               gender * mental_health +
               factor(married) + factor(children) +
               factor(education) + factor(employed) +
               (1 + scale(z_score_age) + factor(married) + factor(children) + 
                factor(education) + factor(employed) | country),
             data = new_data, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))


summary(model)

# Assuming 'model' is a linear mixed-effects model (lmer), and you want to save coefficients to a CSV file

# Extract coefficients from the model summary
coefficients_df <- data.frame(summary(model)$coefficients)

# Write coefficients to a CSV file
write.csv(coefficients_df, "model_coefficients.csv", row.names = TRUE)

```